Google matrix, dynamical attractors and Ulam networks 
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We study the properties of the Google matrix generated by a coarse-grained Perron- Frobenius 
operator of the Chirikov typical map with dissipation. The finite size matrix approximant of this 
operator is constructed by the Ulam method. This method applied to the simple dynamical model 
creates the directed Ulam networks with approximate scale- free scaling and characteristics being 
rather similar to those of the World Wide Web. The simple dynamical attractors play here the role 
of popular web sites with a strong concentration of PageRank. A variation of the Google parameter 
a or other parameters of the dynamical map can drive the PageRank of the Google matrix to a 
delocalized phase with a strange attractor where the Google search becomes inefficient. 

PACS numbers: 05.45.-a, 89.20.Hh, 05. 45. Ac 



I INTRODUCTION 

The World Wide Web (WWW) continues its striking 
expansion going beyond 10^^ web pages. Information 
retrieval from such an enormous database becomes the 
main challenge for WWW users. An efficient solution, 
known as the PageRank Algorithm (PRA) proposed by 
Brin and Page in 1998 [1], forms the basis of the Google 
search engine used by the majority of inter nautes in ev- 
eryday life. The PRA is based on the construction of the 
Google matrix which can be written as (see e.g. [2] for 
details): 



G = aS + (1 - a)^/N . 



(1) 



Here the matrix S is constructed from the adjacency ma- 
trix A of directed network links between N nodes so that 
Sij = Aij/"^^ Akj and the elements of columns with only 
zero elements are replaced by The second term in 

r.h.s. of (1) describes a finite probability 1 — a for WWW 
surfer to jump at random to any node so that the ma- 
trix elements Eij = 1. This term allows to stabilize the 
convergence of PRA introducing a gap between the max- 
imal eigenvalue A = 1 and other eigenvalues A^. Usually 
the Google search uses the value a = 0.85 [2]. By the 
construction Gij = 1 so that the asymmetric matrix 
G has a left eigenvector being a homogeneous constant 
for A = 1. The right eigenvector at A = 1 is the PageR- 
ank vector with positive elements pj and pj = 1. All 
WWW nodes can be ordered by decreasing pj so that the 
PageRank plays a primary role in the ordering of websites 
and information retrieval. The classification of nodes in 
the decreasing order of pj values is used by the Google 
search to classify importance of web nodes. The informa- 
tion retrieval and ordering is based on this classification 
and we also use it in the following. 

It is interesting and important to note that by the con- 
struction the operator G belongs to the class of Perron- 
Frobenius operators [2] . Such type of operators naturally 



appear in the ergodic theory [3] and in the description of 
dynamical systems with Hamiltonian or dissipative dy- 
namics [4, 5]. 

The studies of properties of G are usually done only for 
the PageRank vector which can be find efficiently by the 
PRA due to a relatively small average number of links 
in WWW. At present Google succeeds to operate with 
PageRank vectors of size of the whole WWW being of the 
order of 10^^. It is established that for large WWW sub- 
sets Pj is satisfactory described by a scale-free algebraic 
decay with pj ^ l/j^ where j is the PageRank order- 
ing index and P ^ 0.9 [2, 6]. The studies of PageRank 
properties are now very active in the computer science 
community being presented in a number of interesting 
publications (see e.g. [7, 8, 9] and an overview of the 
field in [10]). 

While the properties of the PageRank are of primary 
importance it is also interesting to analyze the properties 
of the Google matrix G as a whole large matrix. Such an 
analysis can help to establish links between the Google 
matrix and other fields of physics where large matrices 
play an important role. Among such fields we can men- 
tion the Random matrix theory [11] which finds applica- 
tions for a description of spectra in complex many-body 
quantum systems and the Anderson localization which is 
an important physical phenomenon for an electron trans- 
port in disordered systems (see e.g. [12]). A transition 
from localized to delocalized eigenstates also can take 
place in networks of small world type (see [13, 14]). How- 
ever, in the physical systems considered in [11, 12, 13, 14] 
all matrices are Hermitian with real eigenvalues while the 
Perron- Frobenius matrices have generally complex eigen- 
values. 

A first attempt to analyze the properties of right 
eigenvectors ipi (G?/^^ = Xitl^i) and complex eigenvalues 
Xi was done recently in [15]. The Google matrix was 
constructed from a directed network generated by the 
Albert-Barabasi model and the WWW University net- 
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FIG. 1: (Color online) PageRank pj for the Google matrix 
generated by the Chirikov typical map (2) at T = 10, /c = 
0.22, T] = 0.99 (set TIO, top row) and T = 20, /c = 0.3, 
r] = 0.97 (set T20, bottom row) with a = 1, 0.95, 0.85 (left to 
right). The phase space region 0<x<27r;— 7r<p<7ris 
divided on N = 3.6 • 10^ cells; Pj is zero for blue and maximal 
for red. 



works with randomization of links. The Google matrix 
was considered mainly for the value a = 0.85. It was 
shown that at certain conditions a delocalization phase 
emerges for the PageRank and states with complex A. 
In spite of a number of interesting results found in [15] 
a weak feature of models used there is a significant gap 
between A = 1 of PageRank vector and \Xi\ < 0.4 of 
other vectors. We note that according to [15] the Uni- 
versity networks have \Xi\ close to 1 but after random- 
ization of links a large gap emerges in the spectrum of 
A. This gap in |A| was rather large and was not sensi- 
tive to a variation of a in the interval 0.85 < a < 1. 
Hence, the PageRank vector also was not very sensitive 
to a while for real WWW it is know that pj is rather 
sensitive to a due to existence of \Xi\ close to 1 [2, 15]. 
Thus the results obtained in [15] show that even if the 
Google matrix is constructed on the basis of typical mod- 
els of scale-free networks it is quite possible that its spec- 
trum may have a large gap for 0.85 < a < 1 thus being 
rather far from spectral properties of the Google matrices 
of WWW. Therefore it is rather desirable to have other 
simple models which generate a directed network with 
Google matrix properties being close to those of WWW. 

With an aim to have more realistic models we develop 
in this work another approach and construct the Google 
matrix from the Perron- Frobenius operator generated by 
a certain dynamical system. The probability flow in these 
models has rich and nontrivial features of general impor- 
tance like simple and strange attractors with localized 
and delocalized dynamics governed by simple dynamical 
rules. Such objects are generic for nonlinear dissipative 
dynamics and hence can have relevance for actual WWW 
structure. Thus these objects can find some reflections in 



the PageRank properties. The dynamical system is de- 
scribed by the Chirikov typical map [16] with dissipation, 
the properties of this simple model has been analyzed in 
detail in a recent work [17]. We find that the Google 
matrix generated by this dynamical model has many A^ 
close to 1 and the PageRank becomes sensitive to a (see 
Fig. 1). This model captures also other specific prop- 
erties of WWW Google matrices. To construct a net- 
work of nodes from a continuous two-dimensional phase 
space we divide the space of dynamical variables (x^y) 
on N = Nj; X Ny cells (we use = Ny). Then Nc tra- 
jectories are propagated from a cell j on the whole period 
of the dynamical map and the elements Sij are taken to 
be equal to a relative number Ni of trajectories arrived 
at a cell i {Sij = Ni/Nc and Sij = 1). Thus S gives 
a coarse-grained approximation of the Perron-Frobenius 
operator for the dynamical map. The Google matrix G 
of size N is constructed from S according to Eq. (1). We 
use a sufficiently large values of Nc so that the properties 
of G become not sensitive to Nc. 

Such a discrete approximation of the Perron-Frobenius 
operator is known in dynamical systems as the Ulam 
method [18]. Indeed, Ulam conjectured that such 
a matrix approximant correctly describes the Perron- 
Frobenius operator of continuous phase space. For hy- 
perbolic maps the Ulam conjecture was proven in [19]. 
Various types of more generic one-dimensional maps have 
been studied in [20, 21, 22]. Further mathematical re- 
sults have been reported in [23, 24, 25, 26] with exten- 
sions and prove of convergence for hyperbolic maps in 
higher dimensions. However, the studies of more generic 
two-dimensional maps remain rather restricted (see e.g. 
[27]) and non-systematic. In principle the construction 
of directed graphs on the basis of dynamical systems is 
a known mathematical approach (see e.g. [5]) but the 
spectral properties of the Google matrix built on such 
graphs were not studied till now. 

In this paper we show that the Ulam method applied 
to two-dimensional dissipative dynamical maps gener- 
ates a new type of directed networks which we call the 
Ulam networks. We present here numerical and analyti- 
cal studies of certain properties of the Google matrix of 
such networks. 

The paper is organized as follows: in Section II we give 
the description the Chirikov typical map and the way the 
Ulam network is constructed on the basis of this map 
with the corresponding Google matrix, the properties of 
the map and the network are described here; in Section 
III the properties of the eigenvalues and eigenstates of 
the Google matrix are analyzed in detail including the 
delocalization transition for the PageRank, fractal Weyl 
law and the global contraction properties; the summary 
of the results is presented in Section IV. 
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II ULAM NETWORKS OF DYNAMICAL MAPS 
Chirikov typical map 

To construct an Ulam network and a generated by it 
Google matrix we use a dynamical two-dimensional dis- 
sipative map. The dynamical system is described by the 
Chirikov typical map introduced in 1969 for a description 
of continuous chaotic flows [16]: 

yt+i =vyt-^k sm{xt + Ot) , =Xt-\- yt+1 • (2) 

Here the dynamical variables y are taken at integer mo- 
ments of time t. Also x has a meaning of phase variable 
and is a conjugated momentum or action. The phases 
Ot = Ot+T are T random phases periodically repeated 
along time t. We stress that their T values are chosen and 
fixed once and they are not changed during the dynami- 
cal evolution of x^y. We consider the map in the region 
of Fig. 1 (0 < X < 27r, — TT < y < 7t) with the 27r-periodic 
boundary conditions. The parameter < r] < 1 gives the 
global dissipation. The properties of the symplectic map 
at 7^ = 1 have been studied recently in detail [17]. The 
dynamics is globally chaotic for k > kc ^ 2.5/T^/^ and 
the Kolmogorov- Sinai entropy is h ^ 0.29A:^/^ (more de- 
tails about chaotic dynamics and the Kolmogorov- Sinai 
entropy can be found in [3, 4, 28, 29]). 

In this study we use two random sets of phases Ot with 
T = 10 and T = 20. Their values are given in the Ap- 
pendix. We also fixed the dissipation parameter 77 = 0.99 
for T = 10 and 7/ = 0.97 for T = 20. We call these two 
sets of parameters as TIO and T20 sets respectively. The 
majority of data are obtained at A: = 0.22 for the set TIO 
and at A: = 0.3 for the set r20 (see Fig. 1). These are 
two main working points for this work. 

For the set TIO {k = 0.22, r] = 0.99) we have the 
theoretical value of the Kolmogorov- Sinai entropy h = 
0.29^:^/^ = 0.105 for the symplectic map at r] = 1 [17]. 
The actual value at = 1 is determined numerically by 
the computation of the Lyapunov exponent and has a 
value h = 0.0851. For rj = 0.99 we also have the global 
dissipation rate 7c = —Tin 77 = 0.1005 after the map 
period (which is equal to T iterations). The global con- 
traction factor is Tc = r]^ = exp(— 7c) = 0.9043. For a 
weak dissipation the fractal dimension d of the limiting 
set can be approximately estimated in a usual way (see 
e.g. [29]) disd = 2- jc/{Th) = 1.882. 

In a similar way for the set T20 {k = 0.3, r] = 0.97) we 
have the theoretical value h = 0.29k'^^^ = 0.1299, while 
the actual numerical value is = 0.1081. Also here 7c = 
—Tin 77 = 0.609, Fc = 0.5437 and the estimated fractal 
dimension of the limiting set is d = 2 — 7c/ (T/i) = 1.718. 

The bifurcation diagrams for the sets TIO and T20 
are shown in Fig. 2 and Fig. 3 respectively. On large 
time scales we clearly see parameter k regions with simple 
and chaotic attr actors. For a shorter time scales a dis- 
tinction between two regimes becomes less pronounced. 




FIG. 2: Bifurcation diagram showing values of y vs. map 
parameter k for the set TIG of the Chirikov typical map (2). 
The values of y, obtained from 10 trajectories with initial 
random positions in the phase space region, are shown for 
integer moments of time 100 < t/T < 110 (left panel) and 
10^ < t/T < 10^ + 100 (right panel). 




FIG. 3: Same as in Fig. 2 for the set T20. 



This means that during a long time a trajectory moves 
between few simple attractors (which are clearly seen in 
Fig. 1 in the left column) before a final convergence is 
reached. 



Network construction and distribution of links 

The Ulam network for the Chirikov typical map (2) 
is constructed in the following way The whole phase 
space region 27r x 27r is divided into N = x Ny cells 
{Nx = Ny) and Nc trajectories are propagated from each 
given cell j during T map iterations which form the pe- 
riod of the map. After that the elements of matrix Sij 
are computed as Sij = Ni/Nc{j) where Ni is a num- 
ber of trajectories arrived from a cell j to cell i. In this 
way we have by a definition Sij = 1 . Such S gives 
a coarse-grained approximation of the Perron- Frobenius 
operator for the map (2). The Google matrix G of size 
N is constructed from S according to Eq. (1). To con- 
struct Sij we usually use Nc = 10^ but the properties 
of S are not affected by a variation of Nc in the interval 
10^ < Nc < 10^. Since the cell size is very small it is 
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FIG. 4: (Color online) Differential distribution of number of 
nodes with inQoinQ Pin (k) (blue) and outgoing Pout{i^) (red) 
links HI for sets TIO (left) and T20 (right). The straight dashed 



lines give the algebraic fit P(/^) 



with the exponent fi ■ 



1.86, 1.11 (T10,T20) for ingoing and = 1.91, 1.46 (T10,T20) 
outgoing links. Here N = 1.44 • 10^ and P{hv) gives a number 
of nodes at a given integer number of links k for this matrix 
size. Blue point at />: = shows that in the whole matrix 
there is a significant number of nodes with zero ingoing links. 



unimportant in what way Nc trajectories are distributed 
inside the cell. Up to statistical fluctuations, the val- 
ues of Sij remains the same for homogeneous or random 
distribution of Nc trajectories inside a cell. 

Up to = 22500 we used exact diagonalization of G 
to determine all eigenvalues and right eigenvectors i/ji , 
for larger up to = 1.44 • 10^ we used the PRA to 
determine the PageRank vector. The majority of data 
are presented for two typical sets TIO, T20 of parameters 
of the map (2) and the PageRanks for various values of a 
are shown in Fig. 1. For these sets the dynamics has a few 
fixed point attractors but it takes a long time t ~ 10^ to 
reach them. During this time a trajectory visits various 
regions of phase space. 

It is important to note that the discreteness of phase 
space linked to a finite cell size produces an important 
physical effect which is absent in the original continu- 
ous map (2): effectively it introduces an additional noise 
which amplitude a is approximately a ~ 2i\l\fN . This 
becomes especially clear for the symplectic case at 77 = 1 
and at small values of /c at T = 1 (all Qt are the same). 
In this case the map is reduced to the Chirikov standard 
map [28] and the continuous map dynamics is bounded by 
the invariant Kolmogorov-Arnold-Moser (KAM) curves. 
However, the discreteness of phase space allows to jump 
from one cell to another and thus to jump from one curve 
to another. This leads to a diffusion in y direction and 
appearance of a homogeneous ergodic state at A = 1. 
A direct analysis also shows that at any finite cell size 
the operator S has a homogeneous ergodic state with 
A = 1, we also checked this via numerical diagonaliza- 
tion of matrix sizes A^ ~ 20000. This example shows 
that the Ulam conjecture is not valid for quasi-integrable 
symplectic maps in the KAM regime. 

The physical origin of the difference between the con- 
tinuous map and the finite size cell approximation is due 
to introduction of an effective noise term <Jt in r.h.s. of 




FIG. 5: (Color online) Same as in Fig. 3 for the set TIG at 
k = 0.22 (left) (same as Fig. 4 left) and k = 0.6 (right) and 
AT = 3.6 • 10^ The fit gives the exponent /i = 1.87, 1.92 for 
ingoing (blue), outgoing (red) links at A; = 0.22 (left) and \i — 
1.70, 1.83 for ingoing (blue), outgoing (red) links at /c = 0.6 
(right). 



(2) induced by a finite cell size. Due to this noise the 
trajectories diffuse over all region — tt < y < i\ after a 
diffusive time scale Id ~ 11^ j even if the continuous 
map is in the KAM regime with bounded dynamics in 
y. Hence, here a ~ 27i/\/N is an effective amplitude of 
noise introduced by cell discreetness. 

Even if this a-noise leads to a drastic change of dynam- 
ics for quasi-integrable regime its effects are not very im- 
portant in the case of chaotic dynamics where noise gives 
only a small additional variation as compared to strong 
dynamical variations induced by dynamical chaos. With 
such a physical understanding of discreetness effects we 
continue to investigate the properties of the Ulam net- 
works. However, we stress that the <j-noise is local in the 
phase space and hence it is qualitatively different from 
the Google term a which generates stochastic jumps over 
all sites. 

In Figs. 4,5 we show the distributions of ingoing Pin(^) 
and outgoing Pout(^) links n in the Ulam network pre- 
sented by S matrix generated by the map (2) as described 
above. These distributions are satisfactory described by 
a scale- free algebraic decay P ^ 1/ with ji ^ 1.86, 1.11 
for ingoing and 1.91, 1.46 outgoing links at T10,T20 re- 
spectively and a typical number of links per node ^ 10 
(see Fig. 4 and Fig. 5). Such values are compatible with 
the WWW data of scale- free type where ja ~ 2.1, 2.7 for 
ingoing, outgoing links [2, 6]. However, we may also note 
an appearance of certain deviations at large values of k. 
Indeed, for a dynamical system a large number of links 
appears due to exponential stretching of one cell after 
T map iterations that gives a typical number of links 
k ~ exp(/iT). It is possible that during the dynamical 
evolution much larger values of stretching can appear. 
Indeed, the comparison of two cases at /c = 0.22 and 
k = 0.6 for the set TIO in Fig. 5 shows that for larger k 
the scale- free distribution continues to much larger values 
of > 200 while for smaller k the scale-free type decay 
stops around ^ 50. For the set T20 the stretching is 
stronger and the scale-free decay continues up to larger 
values of hz. 
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It is clear that for the Ulam networks discussed here 
one has a rapid exponential decay of links distribution at 
asymptotically large link number k,. However, due to an 
exponential growth of typical ~ exp(/iT) a scale- free 
type decay can be realized up to very large by increas- 
ing T. In these studies we stay at the chosen working 
points where a scale- free decay remains dominant for ma- 
trix sizes of the order of ~ 10^ — 10^. 

Finally we note that the models of the Google ma- 
trix generated by the Ulam networks are most interest- 
ing for dissipative maps. Indeed, by construction the left 
eigenvector of the Google matrix V^^^G = i/j^ dit X = 1 
is a homogeneous vector ip^ = const. As a result for 
symplectic maps the right vector of PageRank pj is also 
homogeneous. Only dissipation term generates an inho- 
mogeneous decay of pj . 



Ill PROPERTIES OF EIGENVALUES AND 
EIGENSTATES 
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FIG. 6: (Color online) PageRank distribution pj for N = 10^, 
9 • lO'^, 3.6 • 10^ and 1.44 • 10^ shown by red, magenta, green 
and blue curves, the dashed straight lines show fits pj ~ 
with p: 0.48 (b), 0.88 (e), 0.60 (f). Dashed lines in panels 
(a),(d) show an exponential Boltzmann decay (see text, lines 
are shifted in j for clarity). Other parameters, including the 
values of a, and panel order are as in Fig. 1. In panels (a),(d) 
the curves at large N become superimposed. Here and below 
logarithms are decimal. 



Delocalization transition for PageRank with a 

The variation of PageRank pj with a is shown in Fig. 1 
for two sets TIG and T20. The distribution pj is plotted 
for each cell of the phase space (x,?/), the numbering 
of cells is done by the integer grid Ux x Uy which has a 
certain correspondence with the index j which numerates 
the values of pj in the decreasing order with j. At a = 1 
the distribution pj is concentrated only on a few local 
spots corresponding to fixed point attractors. Physically 
this happens due to presence of a noise, induced by cell 
discretization, which leads to transitions between various 
fixed points. With the decrease of a the PageRank starts 
to spread over a strange attractor set. The properties of 
strange attractors in dynamical dissipative systems are 
described in [29]. In the map (2) the strange attractor 
appears at larger values of k (namely k > 0.5 for TIO, 
k > 0.34 for T20, see Figs. 2, 3) but a presence of effective 
noise induced by a and 1 — a terms leads to an earlier 
emergence of strange attractor. Below a certain value 
a < ac the PageRank becomes completely delocalized 
over the strange attractor as it is clearly seen in Fig. 1 
for the set TIO. 

The dependence of pj on j is shown in more detail in 
Fig. 6. For a = 1 PageRank shows a rapid drop with 
j that can be fitted by an exponential Boltzmann type 
distribution pj ~ ex.-p{—bjcj/Da) where 6 is a numeri- 
cal constant {b ^ 1.4; 2.1 for T10;T20), 7c = -Tlnr^ is 
the global dissipation rate and D(j = a'^N ^ (^tt)^ is a 
noise diffusion (dashed lines in Fig. 6a, d). Such an ex- 
ponential decay results from the Fokker-Planck descrip- 
tion of map (2) in the presence of a noise term which 
gives diffusive transitions on nearby cells. For a < 1 
random surfer transitions introduced by Google give a 
significant modification of PageRank which shows an al- 



gebraic decay pj ~ 1/j^ with the exponent p dependent 
on a (Fig. 6b,e,f); for the set T20 at a = 0.95 we obtain 
/3 ~ 0.88 being close to the numerical value found for the 
WWW [2]. However, P decreases with the decrease of a 
and for TIO set a delocalization takes place for a = 0.85 
so that Pj spreads homogeneously over the strange attrac- 
tor (see Fig. 1 top right panel and Fig. 6c). For T20 set 
Pj ~ '0i=i(j) remains localized at a = 0.85 so that a PAr- 
ticipation Ratio (PAR) ^ = Ei(IV'i(j)l')V E.dV'iO')!" 
for the PageRank remains finite at large N. We use this 
definion of PAR ^ for all eigenvectors '0i(j). 

Properties of other eigenvectors 

To understand the origin of the delocalization transi- 
tion in a we analyze in Fig. 7 the properties of all eigen- 
values Xi and eigenvectors tjji with their PAR ^. Due 
to a noise activation transitions take place between the 
attractor fixed points leading to states with being ex- 
ponentially close to A = 1 (Fig. 7a). The convergence 
to |A| = 1 is exponential in N for certain states and 
may lead to numerical problems at very large N. How- 
ever, the standard numerical diagonalization methods re- 
mained stable for the values of used in our studies. 

The distribution of A^ in the complex plane is shown in 
Fig. 7c, d: there are A^ approaching A = 1 mainly along 
the real axis but a majority of A^ are distributed inside a 
circle of finite radius around A = 0; this radius decreases 
with the increase of global dissipation from 7c = 0.10 for 
set TIO to 7c = 0.61 for T20. The PAR values for states 
inside the circle have typical values 4 < ^ < 300 shown 
by grayness. The dependence of ^ on 7 = — 21n|A| and 
N shows that the eigenstates inside the circle remain lo- 
calized at large (Fig. 7b). We attribute this to the fact 
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FIG. 7: (Color online) (a) Dependence of gap 1 — |A| on Google 
matrix size N for few eigenstates with |A| most close to 1, set 
TIO, a = 1; (b) dependence of PAR ^ on 7 = -21n|A| for 
N = 2500, 5625, 8100, 10^, 14400 for set TIO, a = 1 (curves 
from top to bottom: red, magenta, green, blue, black); (c) 
complex plane of eigenvalues A for set TIO with their PAR 
^ values shown by grayness (black/blue for minimal ^ ~ 4, 
gray/light magenta for maximal ^ ~ 300; here a = 1, N = 
1.44 • lO'^); (d) same as (c) but for set T20. 




FIG. 8: Probability distribution dW{-f)/d-f for set TIO, a = 
1 at AT = 2.5 • lO^(x), 10^(+), 1.44 • 10^ (dots); W{-f) is 
normalized by the number of states A^^ = 0.55A''^"^^ with 
7 < 6. Inset: dependence of number of states A^ with 7 < 7^ 
on A" for sets TIO (circles, jb = 6) and T20 (triangles, 7^ = 3); 
dashed lines show the fit A^^ = AN"" with A = 0.55, u = 0.85 
and A = 0.97, u = 0.61 respectively. 



that at large N the diffusion due to a noise in presence of 
dissipation leads to spreading only over a finite number 
of cells and thus ^ remains bounded. This ^(7, N) depen- 
dence is different from one obtained in [15] for the Albert- 
Bar abasi model, the comparison with data from WWW 
University networks is less conclusive due to strong fluc- 
tuations from one network to another (see Fig. 4 in [15]): 
an average growth of ^ is visible there even if at A" ~ 10^ 
the values of ^ are comparable with those of Fig. 7b. 
Globally our data of Fig. 7 show that the diffusive modes 
at |Ai| < 1 remain localized on a number of nodes ^ N. 

We also stress an important property of eigenvalues 
and eigenvectors with < |Ai| < 1. In agreement with 
the known theorems [2] our numerical data show that 
for the states with < |Ai| < 1 their are indepen- 
dent of a {Xi are simply rescaled by a factor (1 — a) 
according to [2]). This happens due to a specific prop- 
erty of (1 — a)E/A' term in G, which is constructed from 
a homogeneous vector with rank equal to unity. Right 
eigenvectors are orthogonal to the homogeneous left vec- 
tor and hence (1 — o^) term affects only the PageRank but 
not other eigenvectors. 



Fractal Weyl law for Google matrix 

Another interesting characteristic of G is the density 
distribution dW{'y)/d'y over 7. The data presented in 
Fig. 8 show that its form becomes size independent in the 



limit of large N. At small 7 < 3 the density decreases 
approximately linearly with 7 without any large gap. We 
find rather interesting that the total number of states A'^ 
with finite j < ^ 5 grows algebraically as A"^ = AN^ 
with u < 1 (Fig. 8 inset). We interpret this result on the 
basis of the fractal Weyl law established recently for non- 
unitary matrices with fractal eigenstates (see e.g. [30, 31] 
and Refs. therein). According to this law the exponent u 
isu = d—1 where d is the fractal dimension of the system. 
Approximately we have d—1 ~ 1 — 7c/ {Th) [29, 31] 
that gives u = 0.88,0.72 for the sets T10,T20 with the 
numerical values of 7c, h given above. These values are 
in a good agreement with the fit data u = 0.85,0.61 of 
Fig. 8 inset. The fact that u < 1 implies that almost all 
states have A = in the limit of large A" (in this work we 
do not discuss the properties of these degenerate states 
with large ^ ~ A"). 

It is interesting to note that the fractal Weyl law is 
usually discussed for the open quantum chaos systems 
(see [30, 31] and Refs. therein). There the matrix size 
is inversely proportional to an effective Planck constant 
A" (X l/h. For the Ulam networks generated by dynam- 
ical attractors a cell size in the phase space places the 
role of effective h. This opens interesting parallels be- 
tween quantum chaotic scattering and discrete matrix 
representation of the Perron-Frobenius operators of dy- 
namical systems. 
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PageRank delocalization again 

The dependence of PAR ^ of the PageRank on a and 
N is shown in Fig. 9. It permits to determine the crit- 
ical value ac below which the PageRank becomes delo- 
calized showing ^ growing with N . According to this 
definition we have ^ independent of large N ioi a > ac 
while for a < ac the PAR ^ grows with N . The obtained 
data give ac ~ 0.95, 0.8 for T10,T20. Further investi- 
gations are needed to understand the dependence of ac 
on system parameters. Here we make a conjecture that 
1 — Qfc ^ C^c <^ 1 with a numerical constant C ~ 0.3. 
Indeed, for larger dissipation rate 7c = — Tlnr^ a radius 
of a circle with large density of in the complex plane 
A becomes smaller (see Fig. 7c, d) and thus larger values 
of 1 — a are required to have a significant contribution 
of these excited relaxation modes to the PageRank. Also 
the data of [31] for systems with absorption rate 7c show 
a low density of states at 7 < 7c so that it is natural to 
expect that one should have 1 — o^c ~ 7c to get a signifi- 
cant contribution of delocalized relaxation modes from a 
strange attractor set to the PageRank. It is quite prob- 
able that C depends in addition on system parameters. 
Indeed, even at fixed 7c and a = 0.99 being rather close 
to 1 it is possible to have a transition from localized to 
delocalized PageRank by increasing k in the map (2) (see 
Fig. 9 inset and Fig. 10). This transition in k takes place 
approximately at k ^ 0.55 when fixed point attractors 
merge into a strange attractor (see the bifurcation dia- 
gram in Fig. 2). A peak in ^ around k ^ 0.38 is related 
to birth and disappearance of a strange attractor in a 
narrow interval of /c at /c ^ 0.38. At the same time an 
increase of k from 0.22 to 0.6 practically does not affect 
the link distributions P{i^) changing the value of fi only 
by 10% (see Fig 5). This shows that the correlations in- 
side the directed network generated by the map (2) play 
a very important role. 
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FIG. 9: (Color online) Dependence of PageRank ^ on a for 
set TIG at = 5625 (dotted magenta), 1.44 • 10^ (dotted 
red), 9 • 10^ (dashed red), 6.4 • 10^ (full red) and for T20 at 
iV = 1.44 • 10^ (dotted blue), 9 • 10^ (dashed blue), 6.4 • 10^ 
(full blue). Inset shows dependence of ^ on /c for set TIG at 
a = G.99 with N = 1.44 • IG^ (dotted red), 9 • IG^ (dashed 
red), 3.6 • IG^ (full red). 




FIG. IG: (Color online) Same as Fig. 1 for the set TIG at 
a = G.99, = 3.6 • IG^ at /c = G.22 (left) and k = G.6 (right); 
PAR ^ are the same as in the inset of Fig. 9. 



Global contraction 

As discussed above a nontrivial decay of the PageRank 
in our Ulam network appears due to a dissipative na- 
ture of the map (2). Indeed, since 7^ < 1 there is a global 
contraction of the phase space area by a factor Fc = rf 
after T iterations of the map (after its period). Such a 
property is very natural for the continuous map but it is 
more difficult to see its signature from the matrix form 
of the Perron-Frobenius operator after the introduction 
of discreteness of the phase space. 

Nevertheless this contraction can be extracted from the 
matrix G taken at a = 1. To extract it we apply G with 
a = 1 to a homogeneous vector p^-^^ = 1/7V getting the 
new vector p^^^ = Gp^^^ and count the number of nodes 
A^r where p^^^ > qjN and < g < 1 is some positive 



number characterizing the level of the distribution. Then 
the contraction of the network is defined as a fraction of 
such states: F TVp/A^- 

The result of computation of the contraction factor for 
the Ulam network of map (2) for the sets TIO, T20 is 
shown in Fig. 11. The network contraction parameter F 
is independent of g in a large interval 10 ~^ < q < 0.1 and 
it converges to the contraction value Fc of a continuous 
map in the limit of large matrix size N. 

We think that the Google matrix of WWW networks 
can be also characterized by a global contraction factor 
and it would be interesting to study its properties in more 
detail. However, this remains a task for future studies. 
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FIG. 11: (Color online) Dependence of the network contrac- 
tion factor r on the level q of probability distribution over the 
network nodes (see text). Left panel shows data for the set 
TIO at k — 0.22, right panel shows data for the set T20 at 
k — 0.3 for the Ulam network of map (2). The size of the net- 
work is = 10^, 4- 10^, 16 • 10^ (curves from top to bottom at 
q — 0.01). The dashed curves show the contraction Fc = r]^ 
of the continuous map (2) corresponding to the network with 
iV = oo. 



IV SUMMARY 

In summary, we demonstrated that the Perron- 
Frobenius operator built from a simple dissipative map 
with dynamical attractors generates a scale-free directed 
network with properties being rather similar to the 
WWW. The networks and their Google matrices are 
obtained on the basis of the Ulam method for coarse- 
graining of the Perron-Frobenius operator and thus can 
be viewed as the Ulam networks or Ulam graphs. The 
Google matrix of such Ulam networks reproduces many 
properties of real networks with an algebraic decay of 
the PageRank and quasi-degeneracy of eigenvalues near 
unity for the Google parameter a = 1. In this formula- 
tion the popular websites correspond to dynamical fixed 
point attractors which help to generate global scale-free 
properties of the network. The PageRank of the system 
becomes delocalized for a smaller than a certain criti- 
cal value, such a delocalization is linked to emergence 
of a strange attr actor. Even for a very close to unity 
a moderate change of system parameters can drive the 
system to a strange attractor regime with a complete de- 
localization of the PageRank making the Google search 
inefficient. In view of a great importance of the Google 
search for WWW [2, 6] and its new emerging applica- 
tions [32] it may be rather useful to study in more detail 
the properties of the Google matrix generated by simple 
dynamical maps. 

Of course, it is quite possible that at the present state 
the Google matrix of WWW is more stable in respect 
to variation of a (indications for that can be found e.g. 
in [8, 9, 10]). However, WWW evolves with time and 
may become more sensitive to changes of a. Also the 
Google search can be applied to a large variety of other 
important networks (see e.g. [10, 32]) which may be more 
sensitive to various parameter variations. It is quite pos- 
sible that the Ulam networks discussed here only par- 



tially simulate the properties of the WWW. However, the 
Ulam networks are easy to generate and at the same time 
they show a large variety of rich interesting properties. 
The parallels between the Ulam networks and the actual 
WWW can be instructive for deeper understanding of 
both. Therefore, we think that their further studies will 
give us better understanding of the Google matrix prop- 
erties. The studies of the Ulam networks will also lead 
to a better understanding of intricate spectral properties 
of the Perron-Frobenius operators. The application of 
the thermodynamical formalism [33, 34] to the spectra 
of such operators can help to understand their properties 
in a better way. 
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APPENDIX 

The Chirikov typical map (2) is studied here for the 
following random phases ^t/27r for the set TIG : 
0.562579, 0.279666, 0.864585, 0.654365, 0.821395, 
0.981145, 0.478149, 0.834115, 0.180307, 0.15902 
and for the set T20: 
0.415733267627, 0.310795551489, 
0.924301928270, 
0.647524548037, 
0.370553510280, 
0.834758628241, 
0.828693568678, 
0.278538720979, 0.653773338142. 

The numbers are ordered in the serpentine order for 

i = i,2,...r. 

After each T iterations the values of y are reduced 
inside the interval (— tt, tt) corresponding to the periodic 
boundary conditions. 



0.749488203411, 
0.118768635110, 
0.952312529146, 
0.814808044380, 
0.702057578688, 



0.632094907846, 
0.635937571045, 
0.651928927275, 
0.810837257644, 
0.993694010264, 
0.855421638697, 



